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April 17, 2013 

Abstract Granular segregation is an important mech- 
anism for industrial processes aiming at mixing grains. 
Additionally, it plays a pivotal role in determining the 
kinematics of geophysical flows. Because of segregation, 
the grainsize distribution varies in space and time. Ad- 
ditional complications arise from the presence of com- 
minution, where new particles are created, enhancing 
segregation. This has a feedback on the comminution 
process, as particles change their local neighbourhood. 
Simultaneously, particles are generally undergoing remix- 
ing, further complicating the segregation and comminu- 
tion processes. The interaction between these mecha- 
nisms is explored using a cellular automaton with three 
rules: one for each of segregation, comminution and 
mixing. The interplay between these rules creates com- 
plex patterns, as seen in segregating systems, and depth 
dependent grading curves, which have been observed in 
avalanche runout. At every depth, log- normal grading 
curves are produced at steady state, as measured ex- 
perimentally in avalanche and debris flow deposits. 
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1 Introduction 

Flowing granular material exhibits complex behaviour, 
including many phenomena that cannot be described 
in the context of a conventional fluid. For example size 
segregation, comminution and agglomeration all have 
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Fig. 1 Top: A mill stone — the grinding of grain using a mill 
stone is one of the oldest industrial problems in human his- 
tory, yet still mathematically unsolved. Particles of grain are 
crushed to a fine powder by very large deformation shearing 
at high normal stress. The fine powder segregates out of the 
shear zone into cavities built into the mill stone, and then 
under the action of centripetal forces migrates out to a col- 
lection bin. Bottom: A long run-out landslide, where the ratio 
of L/H can be up to 10. L and H are the change in horizontal 
position and height respectively of the centre of mass of an 
avalanche during run-out. Large values of L/H are possible 
indicators of lubrication by a layer of very small particles at 
the base of the flow, which have been created as a result of 
comminution and percolated downwards through segregation. 



no analogue in traditional fluids. To describe these phe- 
nomena the grainsize distribution has o be involved as 
a dynamic property. 

The dynamics of granular material are important in 
many natural processes, such as debris flows, landslides, 
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rockfalls and shear banding. Industrial processing also 
requires many granular flows, such as rotating or tum- 
bling mills, chute flows and hopper filling or discharge. 
In these types of flow processes, breakage of particles 
can be important. As yet, there are no temporo-spatial 
continuum models for measuring changes in the grain- 
size distribution for such open systems where particles 
can advect in space. 

This is especially important in two poorly under- 
stood systems, shown in Figure [T] The first is ancient 
grain milling, where combined normal and shear stresses 
crush wheat grains, dynamically sieving the resultant 
mix. The second is long run out avalanches, where it is 
unclear why these incredibly destructive natural phe- 
nomena can travel enormous distances, up to 10 times 
their vertical fall [HE]. 

With regards to natural flows, it has been under- 
stood for some time that there is a need to model spatial 
and temporal variability in the grainsize distribution of 
flowing material to be able to implement appropriate 
rheological models [3]. 

Here we tackle this problem using a cellular automa- 
ton [4] with three distinct rules of operation: segrega- 
tion, remixing and comminution. As has been shown 
previously [5] for bidisperse systems, we can describe 
segregation in terms of the swapping of cells in a cellu- 
lar automaton. Comminution rules have also been de- 
veloped [6j[7], but these are limited to closed systems. 
Here we will present rules which can apply in open sys- 
tems to arbitrarily polydisperse materials. 

As in [8 we denote the grainsize, s, as an inter- 
nal coordinate of the system such that every point in 
space has a grainsize distribution. We then describe this 
continuous grainsize distribution (j)(s) of the system in 
terms of the solid fraction <P(s) of particles between 
grainsizes s a and s^ as 



<P[s a < s < s h ] = / 0(s') ds'. 

J s n 



(i) 



Conservation of mass at a point in space r = {#, y, z} 
can then be expressed as 



where material can move spatially. The right hand side 
of the same Equation represents mechanisms tradition- 
ally treated as closed systems, such as agglomeration, 
crushing and abrasion. Each of these systems — open 
and closed — has been the subject of much study, but 
the coupling of such processes using a continuum de- 
scription has yet to be achieved. 

2 Cellular automata and continua 

To model these systems we use a cellular automaton 
in two spatial dimensions. It is a series of cells on a 
2D regular cartesian lattice of size N x by N z , with di- 
rections x and z, and cell spacing Ax and Az in the 
respective directions. The z direction is perpendicular 
to the shear direction, such that for inclined plane flow 
it points normal to the slope. The x direction represents 
a micro scale internal coordinate. 

We allow the grainsize distribution to be a function 
of this internal spatial coordinate such that s = s(x). 
By discretising the grainsize distribution cj)(s) into N x 
monodisperse components of equal volume, they can be 
arranged in the x direction such that summing over this 
coordinate would recover the full grainsize distribution. 
We consider the local neighbourhood of a particular 
particle as those that are adjacent in the x direction. 
The x direction now contains more information than 
the grainsize distribution alone, as the local orientation 
of particles is preserved, below the resolution of the 
analogous continuum scale. 

We number the cells from the bottom-left corner of 
the grid so that position on the grid can be expressed 
using the pair {i, j}, where i and j indicate the number 
of cells across in the respective z and x directions. In 
all cases the system is considered to be periodic in the 
j direction. 

Each cell contains a single number, Sij, which dic- 
tates the grainsize of the particles in the representative 
volume element defined by the cell {i,j}. 

We can define a discretised grainsize distribution <pi 
at any height i as a histogram of the number of cells 
within a discrete grainsize fraction with centre s a and 
width As in all N x neighbours taken in the j direction: 



dt 



V-(0u) = /i + -/i-, 



(2) 



where u(r, s,i) = {u : v,w} is the material velocity, 
/i + (r, s,i) is the birth rate, describing the creation of 
new particles of grainsize s at time t, and ft~(r, s,t) 
is the death rate, at which particles of grainsize s are 
destroyed. 

The second term in Equation [2] describes the ad- 
vection of mass, such as characterises open systems, 



l{Sa) ~ N x As ^ 1 
x k=i I 



1 if Sa -^< Si|fc < Sa + ^, 
otherwise. 



(3) 



We also define the local average grainsize over the 
nearest neighbours in the j-direction as 



&i,j — \Si,j — l ~r Sij+i)/ z, 



(4) 
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Fig. 2 The cushioning effect and nearest neighbour rule. Left: 
Large particles are cushioned such that they will not break 
because of an abundance of small particles. Right: Small parti- 
cles do not carry a significant amount of load as they are free 
to move in interstitial pore spaces. Top: Initial conditions. 
Bottom: Result after one iteration of the cellular automata 
rule, where on the left the small particles in a cell (orange) 
become even smaller (green), and on the right the large par- 
ticles in a cell (blue) become smaller. 



Following the formulation in [10] , for the case where 
particles are created only by the fragmentation of larger 
particles, we can express the death rate as 



h-(s,t) = b(s)cj)(s,t), 



(5) 



where b(s) is some specific breakage rate which governs 
the frequency at which particles of grainsize s break 
into smaller fragments. The birth rate is then the sum 
of breakages into size s over all particles larger than s, 
which can be expressed as 



/•OO 

h + (s,t) = / b(s')P(s\s')(l>(s',t) ds', 

J S 



(6) 



where P(s\s') is a probability density function which 
dictates the probability of creating grainsize s from 
crushing a particle of grainsize s'. We can then express 
conservation of mass as 



dt 



b(s')P(s\s')<l)(s',t) ds' -b(s)(j)(s,t). (7) 



In a discrete sense, such as that defined in the cel- 
lular automaton, we can rewrite this equation as the 
conservation of a grainsize fraction with centre s a and 
width As, over a time step At as 



and similarly for the i-direction. For the cellular au- 
tomata rules defined below, we need to define two new 
operators: s a O 55 represents the swapping of values 
s a and Sfo between their respective cells, and s a =>• s^ 
represents a change in grainsize in the cell containing 
s a to the new value of s&. 

In the limit of infinitesimal cell size, the cellular au- 
tomaton may be regarded as a first order partial differ- 
ential solver [11 J. In fact, there are many ways to recover 
continuous results from such a cellular automaton [12] . 



3 Closed systems 

We begin by considering closed systems, which are those 
in which material does not advect in space such that 
u = 0. There are many processes that have previously 
been represented as closed systems, such as comminu- 
tion, agglomeration and abrasion. We here look only 
at comminution, where particles are crushed to form 
fragments of smaller sizes. 



A(f)j(s a ) 

At 



N s 

^2 (bi(sk)P(s a \sk)(l>i(sk)As 



k=a+As 



bi(s a )<f>i(8 a ), 



(8) 



where N s is the total number of evenly spaced bins 
of size As. These equations have been considered many 
times before [T3j[l4l[T5, 16 , and solutions have been pro- 
posed for many mechanisms of comminution, such as 
grinding, cleavage and abrasion. However, breakage mech- 
anisms are normally assumed with a priori knowledge 
of power law distributions [9|fT3] . In fact, in most mod- 
els, either the breakage rate 6, the fragment probability 
distribution P, or both, are generally assumed to be 
power law in nature from the outset [T3] . 

Another method to model the problem has been 
proposed in various forms, and uses simple geometric 
analogies in a cellular automaton [6,7 where power law 
patterns are found, not imposed, by assuming that par- 
ticles with neighbours of the same size are likely candi- 
dates to crush. 

We can unify these two approaches, of macroscopic 
grainsize distribution changes, and microscopic near- 
est neighbourhood behaviour, by including the grain- 
size distribution in a cellular automaton. 
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Fig. 3 Evolving cumulative grainsize distributions due to 
comminution. All simulations have N z = 1, N x = 10 6 and 
kb = 1. Three initial conditions are considered, with oti = -5, 
-2 and 1. Each initial condition is crushed at three values of 
(3 = 0.05, 0.2 and 0.8. Within each sub-plot, each line repre- 
sents a different time t = 0, 0.05, 1, 2 and 10, from bottom 
right to top left. 



3.1 The crushing mechanism 

When a large particle is surrounded by small particles, 
it is cushioned from fracture by having many points of 
contact with neighbours, leading to a fairly isotropic 
loading state, as shown on the left of Figure [2] Also, 
when a small particle is surrounded by large particles, 
it does not carry significant load, as it is either able to 
fit in the pore spaces if sufficiently small, or is highly 
mobile, as on the right of Figure [2] Because of this, 
we consider fracture of particles only when they are 
surrounded by particles of a similar size. We then have 
a condition for fracture such that for a particle of size 
Si j and a local neighbourhood of particles of average 
grainsize Sij, 



(0, 1) X 8 id 



if 



i\</3sij 



(9) 



This rule states that all particles in the representa- 
tive volume defined at Sij reduce in size by a randomly 
chosen factor between and 1 if it is within /3sij of 
the local mean grainsize. The non-dimensional factor (3 
determines how similar particles must be to their neigh- 
bours before crushing can occur. We have also included 
a factor of Sj on the right hand side of the inequal- 
ity to make small particles harder to crush, recognising 
that smaller particles have a higher crushing stress than 
larger ones [17] . This crushing event has some frequency 
which is proportional to the shear strain rate, allowing 
us to define the breakage rate as 




10° 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 
Normalised grainsize, s/s m 

Fig. 4 Cycles of breakage towards an attractor. Top: Cumu- 
lative grainsize distributions at different numbers of cycles 
of loading. Each line represents the grading after a certain 
number of crushing-remixing cycles. From bottom right to 
top left, they are 0, 1, 10, 50, 100 and 500 cycles. The system 
approaches a = 3, or a horizontal line in this plot. Bottom: 
The same data as in the top plot, now shown as number of 
particles A greater than a certain size s. Insets show zoomed 
areas of original plot. 



where fc& is a non-dimensional fitting parameter and H 
is the Heaviside function. In ^ we have also defined 
the fragment size distribution as 



P(s a \s hJ )=As/s K 



j> 



(ii) 



where As is the size of the grainsize bin in the s di- 
rection containing fragment size s a . The simulation oc- 
curs over N x cells, spaced Ax apart. Initial conditions 
are generated by sampling N x times from F = Ax to 
F = 1 linearly along the inverse power law distributions 
defined by 



pi/(3-ai) m 



(12) 



bij = h\^i\n(/3s 



1,3 



D i,3 I 



(10) 



where F(s) = J s (j)(s') ds' is the cumulative grainsize 
distribution function and 3 — c^ is the power law gra- 
dient. By time marching with a sufficiently small time 
step, such that hi At < 1, we implement the frequency 
of breakage in a probabilistic manner. 

Figure [3] shows three different initial conditions, and 
their progression to a steady state grainsize distribu- 
tion. Each initial condition is considered for three val- 
ues of /?. For each case, a new distribution is reached 
after a single timstep of At = 0.05. As time progresses, 
(3 controls most of the behaviour of the process. For 
small values of /?, a steady state is reached where the 
grainsize distribution does not change appreciably af- 
ter time t = 5, resulting in a power law dimension of 
a = 1.99 ±0.01. For large /?, all of the largest parti- 
cles are continually crushed, eventually resulting in a 
system where a approaches 3. 
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L, Length of interrogation window 

Fig. 5 Fractal dimension a against interrogation window 
length L. a value is measured as a best fit of the cumula- 
tive grainsize distribution from s = 10 -4 to 10 _1 , restricted 
to the cells between j = a±L/2. Each line represents a differ- 
ent value of a = 0, N x /4, N x /2 and SN x /4. For this simulation 
N x = 10 6 , (3 = 0.2, k b = 1 and a { = -2. 



The cellular automata predicts the same final grain- 
size distribution largely independent of the initial grad- 
ing, with some minor effects due to the initial concen- 
tration of very large particles. Such a power law grad- 
ing of the grainsize distribution has been measured in 
other cellular automata [6,7J, discrete element simula- 
tions [18] and experiment [2]. 

Generally fractal dimensions are measured in fault 
gauges, confined comminution tests and rock avalanches 
in the range of a = 2 - 3 [T9ll2Qll2T] . The limiting value 
of a = 2 for most cases of our model can be explained 
using the cellular automaton developed in [20] . where 
every particle in the systems has the same probability of 
crushing, given some additional geometric constraints. 
If this probability is exactly 0.5, the system develops 
a fractal dimension of a = 2. If the probability is 1, 
the system reaches a = 3, which represents a system 
with large strain, where nearest neighbours change over 
the crushing period [22]. This distribution in fact cor- 
responds to a random appollonian packing, as in [23] . 

In Figure [5] a is measured many times in a single 
simulation. A number of cells, L, centred at a, are cho- 
sen, and a best fit of a is measured for the cumulative 
grainsize distribution. With increasing zoom into the 
system, the measured power law does not change slope 
significantly. This is an indication of a fractal distribu- 
tion [6 . 

The idea of such final power law grainsize distribu- 
tions has been applied in [24 to develop a thermody- 
namically consistent theory for comminution processes 
in granular materials in closed systems. Extension of 
this or other theories to open systems where particles 
can advect has not yet been achieved. 
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Fig. 6 Bidisperse cellular automaton. Top: Schematic repre- 
senting bidisperse segregation in 2D flow down an inclined 
plane. Bottom: The complimentary ID cellular automaton, 
where large particles and small particles swap over time. 



3.2 Comparison with continua 

To compare the cellular automata rules with a contin- 
uum, we consider the evolution of a large number of 
cells simultaneously, and find averaged properties that 
represent the continuum scale. We express the breakage 
rate, 6^, of a single grainsize fraction s a covering sizes 
over a range of As as 



hi 






N x 



^$>(/fc 



M 



^J 



^,J 



I). 



(13) 



To solve this system globally, we need to sum over 
j, which represents local information about the nearest 
neighbours. If we were to represent this in a continuum 
sense, this local information is smaller than the contin- 
uum scale, and so a new length scale, (", must be intro- 
duced. In the cellular automaton, the length scale which 
controls this behaviour is that of the nearest neighbour 
zone (here set arbitrarily to Ax) . By increasing the size 
of the neighbourhood over which we find the local av- 
erage grainsize, the system would converge towards a 
different state, where physical proximity of neighbours 
is not considered. 

How to introduce such a length scale in a continuum 
theory is an open question. In the cellular automaton, 
we can consider changing this length scale by doing 
cycles of crushing, as will be shown below. 
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Fig. 7 The segregation mechanism as a function of gransize s. 
Left: Bidisperse rule used previously in [5] Right: Polydisperse 
rule used here, where segregation frequency is a function of 
distance to the local mean size. 



3.3 Cycles of crushing: Towards open systems 

We can extend this simulation to a quasi-open system 
by considering a rearrangement of particles within the 
cellular automaton, but without true advection. 

After the system has reached steady state, and no 
further significant crushing will occur, we shuffle the 
system, relocating all of the nearest neighbours, and 
resume crushing until a new steady state is reached. We 
can continue this cycle until an ultimate steady state is 
reached. 

We begin with the simulation shown in the centre of 
Figure |3j which has an initial grading defined by on = 
—2, and after one full crushing iteration process has 
reached a = 2. As progressively more iterations occur, 
a region of higher slope develops at larger grainsizes, 
and this propagates to lower grainsizes with increasing 
iterations, as shown at the top of Figure [4] This trend is 
towards a = 3, however this coincides with a horizontal 
line in this plot, and cannot be seen clearly. To visualise 
this more clearly, the bottom plot in Figure [4] shows the 
same data, but plotted as the number of particles A 
greater than a certain size s, such that the slope of the 
graph is a. We define the number of particles per cell as 
being inversely proportional to the size cubed. The final 
state has been shuffled and crushed 500 times, tending 
towards an ultimate grading with this new power law 
gradient of af =3. 

This effect of cycles of crushing towards af = 3 has 
been observed experimentally [25], numerically with a 
crushable discrete element method [26] and predicted 
analytically as the maximum entropy path towards the 
least efficient packing of the system [23] . It is remark- 
able that such a simple, one dimensional system as this 
can replicate the packing involved in such a complex 
system. 



4 Open systems 

In order to model open systems, where particles can 
advect between points in space, we need rules for ad- 
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Fig. 8 The time evolution the average grainsize s of a bidis- 
perse simulation with varying N x subject to segregation only. 
All cases have k s = 1 and N z = 100. Clockwise from Top Left: 
N x = 1, 10, 100 and 1000. For the case of N x = 1 we find the 
local average grainsize in the z direction. Bottom right: Black 
lines indicate positions of concentration shocks from solution 
of analogous continuum equation. 



vection. To be a truly physical model, we need to satisfy 
the conservation equation for grainsize (j) such that 



dt 



V • (0u) = 0. 



(14) 



We set up our rules for the cellular automaton such 
that this conservation law is held for some velocity u. 
The first mechanism we consider is due to segregation. 
Towards this end, we begin with the simplest case of 
segregation and describe bimixtures, where we present 
a model similar to that proposed in [5]. 



4.1 The segregation mechanism 

As particles flow they collide, creating new void spaces 
which are preferentially filled by smaller particles mov- 
ing, as in Figure [6] The rate of creation of void spaces 
is governed by the shear strain rate, 7. 

The simplest description of this system in terms of 
grainsize is shown on the left of Figure [7] for a bimixture 
of sizes s m and %, where the swapping frequency / is 
defined such that small particles always move down, 
and large particles move up. 

We can extend this model to describe polydisperse 
materials by using the formulation developed in |5], 
where it was shown through energy considerations that 
if a particle is larger than the local average, it has some 
probability of moving up, and conversely if it is smaller 
than the average it will move down. Increasing distance 
in the s direction from the average will increase the like- 
lihood of swapping linearly. This is shown on the right 
of Figure [f| 
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Fig. 9 Time evolution for bidisperse shear flows subject to 
segregation. For all cases k s = 1. Left to right: The system 
is initially filled with 20%, 50% and 80% small particles re- 
spectively. Top row: Simple shear, where 7 = 1. Bottom row: 
Inclined plane flow shear condition, where 7 = VI — z. 



We facilitate this movement by swapping this grain- 
size with that either above, s*+i,j or below, s*-i,j de- 
pending on whether it is smaller or larger than the 
average size Sij. We define the rate of swapping as 



/ = &s 1 7* In- 
direction: 



n 



'1,3 



1), with the sign determining the 



if 



/ > 0, s id O s i+1J 
f < 0, Si j <&Si- lt j. 



where k s is a non-dimensional parameter controlling the 
rate of segregation. We iterate in two half time steps, 
alternately applying this rule firstly to all odd rows, 
and then all even rows, so that particles are inhibited 
from moving very large distances in a single time step. 
Additionally, we only allow swapping upwards if the 
particle is larger than the one above it, or smaller than 
the one below it if moving downwards. 
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Fig. 10 Polydisperse segregation under simple shear. For all 
cases k s = 1 and 7=1. Top to Bottom: Each row represents a 
single simulation with initial condition defined by a^ = -2, 
and 2 respectively. Left: Initial cumulative grainsize distribu- 
tion at three different heights. The solid red, dashed green and 
dash-dotted black lines represents z = 0.9, 0.5 and 0.1 respec- 
tively. Middle: Plot of the average grainsize s over height and 
time. Right: Final cumulative grainsize distributions, plotted 
in the same manner as the initial grainsize distributions. 



rate profile predicted in [8] for the case of inclined plane 
flow: 7 = y/(l — z)/s. In each case, complete segrega- 
tion is observed, where every large particle lies above 
every small particle. The non-uniform shear strain rate 
in the bottom case causes non-uniform transient be- 
haviour towards a steady grading that is the same as 
the top case. According to the formulation presented 
here, the final grading is independent of the loading 
condition. 

As shown in [5 , this model represents a very simple 
analogy of the analytic works done by [28| f27] to model 
bidisperse segregation, with a very simple extension to 
polydisperse systems, as shown below. 



4.2 Bidisperse segregation 

To model a simple bidisperse material, we take a single 
column of a bimixture (N x = 1), of equal proportions of 
sizes s m and sm, randomly allocated to cells, and allow 
it to segregate under simple shear with 7(2) = 1 and 
k s = 1. The result is shown in the top left of Figure 
|8] We can then run the simulation with N x = 5 and 
average over the x-direction. The result of this is shown 
in the top right of Figure [8] We can do this repeatedly, 
to get an increasingly resolved image of the process, as 
shown in the bottom row of Figure [8] for N x = 50 and 
1000. With increasing resolution, this converges on the 
analytic solution presented in [8] and [27] . 

Figure [9] shows the average grainsize at each height 
over time for two different shear regimes, each for 3 dif- 
ferent initial concentrations of small particles, s m . The 
top row depicts simple shear, as in Figure [8j while the 
bottom row uses a simplified version of the shear strain 



4.3 Polydisperse segregation 

We can create a polydisperse sample by generating ini- 
tial conditions in the same way as previously described 
for the breakage cellular automaton, using Equation [T2| 
The segregation patterns produced for a range of ini- 
tial conditions at constant segregation rate k s = 1 and 



with 7 = 1 are shown in Figure 10 Since this is now a 
polydisperse sample, we can calculate the grainsize dis- 
tribution (j)(s). On the left hand side of Figure 10 are the 



initial cumulative grainsize distributions, which are ho- 
mogeneous. During the simulation, segregation occurs, 
creating a non-homogeneous steady state condition af- 
ter some time. These grainsize distributions, which now 
vary with height, are shown on the right hand side of 
the same Figure. 

Averaging over all cells at a given height z, we can 
express the mean segregative velocity Ui of a single 
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Fig. 11 The mixing mechanism. Ten cells, initially segre- 
gated with all large particles (blue) above small particles (yel- 
low), subjected to the mixing mechanism only. Over time, the 
system reaches a disordered state. 



grainsize fraction centred at s a from all cells at height 
i in a time At as 




Time, t 

Fig. 12 The time evolution the average grainsize sofa bidis- 
perse simulation with varying N x subject to mixing only. All 
cases have D — 0.01 and N z = 100. Clockwise from Top Left: 
N x = 1, 5, 50 and 1000. Bottom right: Solid lines indicate 
contours from solution of analogous continuum equation. 



This model, together with the rule for remixing, 
which will be shown next, represents the simplest de- 
scription of the analytic description presented in [8] . 



N x 



Ui(s a ) 



N n 



^2fij(Sa) 






N x 



(15) 



We have included the local average grainsize sij in 
the formulation so that we know which particles are 
locally small or large. As particles are being swapped 
between heights — between representative volume ele- 
ments at the continuum scale — we only require a single 
average grainsize per height, and can in this case freely 
extend the neighbourhood domain over which we find 
the average grainsize s to include every cell at height 
j, labelling it now Si = 1/N X V • Sij- In this case, the 
mean velocity can be expressed as 



4.4 Remixing 

In nature we rarely see such perfect segregation as that 
pictured above. This is due to the random fluctuation of 
particles as the flow propagates down slope. As has been 
done before analytically [30], we can capture this effect 
by introducing remixing into the flow. For the simplest 
case, we allow particles to swap randomly either up or 
down with some frequency, given by a constant D/Az 2 . 
At this stage we let this probability be independent of 
the shear strain rate 7, although a strong dependency 
has been observed [31] in experiments. With frequency 
of swapping controlled by the diffusivity, D, 



Ui(s a ) = hh 



" l Si 



1 



(16) 



Compare this with the analytic description of the 
segregation velocity with no diffusion as predicted by 
[8] for a continuum with internal grainsize coordinate 



s, 



. N . . . a cos 6 ( s 

u(* = 7" = 

c \s 



(I- 1 ). 



(17) 



where c is a fitting parameter, g is the acceleration due 
to gravity, and 6 is the angle of the plane down which 
flow is occurring. As shown in [29], cellular automata 
can successfully be used as a coarse finite differencing 
method to model systems such as these without resort- 
ing to complicated flux limited finite difference schemes, 
as would otherwise be necessary. 



Flip a coin, if heads: Si <^> si-i 
if tails: s^ O s i+ i 

An example of the mixing rule acting on a single 
column of cells over time is shown in Figure [TT] Ini- 
tially, the system is perfectly segregated, but over time 
the systems becomes randomised due to the presence of 
remixing. The characteristic time for mixing to occur is 
the inverse of the diffusivity D/H 2 . 

We are describing a system of cells undergoing Brow- 
nian motion, whereby particles move by the application 
of random forces over time scales that are short rela- 
tive to the motion of the particle. When considered over 
long time scales and large numbers of particles, this is 
analogous to Fickean diffusion [31] . Many other cellular 
automata exist to model pure diffusion [32] . 

As in the case of segregation, this process can be 
averaged over the x direction to describe the evolution 
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Fig. 13 Coupled comminution and mixing. Varying values 
of D/(k b H 2 ). For all cases, N z = 21, N x = 1000. Plot shows 
final value of best fit to power law part of cumulative grainsize 
distribution at t = 5, 50 and 500, corresponding to the solid, 
dashed and dash-dotted lines respectively. 



of the average grainsize at any height over time. By 
increasing the number of cells in the x direction, we 
can increase the smoothness of our solution. Figure 12 
shows the same system as Figure [TT] initially segre- 
gated, that mixes over time to create a homogeneous 
system, but now for increasing numbers of cells in the 
x direction. 

This diffusive behaviour can be described at the con- 
tinuous limit using Fick's first law of diffusion, 







o 



Time, t 



Fig. 14 Bidisperse segregation under simple shear with diffu- 
sion. Left to Right: Increasing diffusivity D = 0.002, 0.01 and 
0.05 with constant segregation coefficient k s = 1. Increasing 
the mixing coefficient smooths out concentration shocks in 
the spatial direction, giving more physically representative 
solutions. 



tension of the cycles of crushing pictured in Figure |4j 
but now with true advection. In this case, as in Figure 
[4j we expect that after a short time relative to the dif- 
fusive time, the system will reach af = 2, as significant 
mixing has not yet occurred. At longer times, the sys- 
tem will approach af = 3, and its final grading. This 
effect is captured in Figure [l3j where the diffusive time 
is controlled by the ratio D/(ki)H 2 ), and the system is 
constrained at af = 2.91. 

For vanishingly small diffusivities, the system will 
still approach af = 3, but only after very long periods 
of comminution. Conversely, at very large diffusivities, 
the system passes af = 2 very rapidly, and approaches 
af = 3 in a relatively short time. 



dt 



D 



d 2 ^ 
dt 2 ' 



5 Coupled problems 



We now have three distinct processes which can be de- 
scribed simultaneously in a single simulation. These 
have all been shown above with their analogous con- 
tinuum description, yet not in all cases could a direct 
link be shown. For the case of comminution, an internal 
length scale governing the spatial distribution of grain- 
size over a sub-continuum length scale was required. 

As all of the mechanisms previously described have 
been created in the same framework, we can simply run 
a cellular automaton which includes multiple phenom- 
ena at the same time. As will be shown in this Section, 
we can investigate the interactions between the mech- 
anisms by varying the parameters which control their 
effects. 



5.1 Comminution and mixing 

We begin with a cellular automata that includes both 
comminution and mixing. This system represents an ex- 



(18) 5.2 Segregation and mixing 



In flows of polydisperse granular materials where com- 
minution does not occur, we can model the evolution of 
the grainsize distribution as being comprised of segrega- 
tive and diffusive remixing components. This occurs in 
many industrial mixing processes, and may be sufficient 
to model levee formation and runout characteristics in 
landslides. At higher speeds remixing increases, sup- 
pressing segregation, while at low speeds segregation 
can play a dominant role in the flow behaviour. 

The effect of coupling mixing and segregation in a 
bimixture can be seen in Figure [l4j where increasing 
diffusivity D smooths out the concentration shock be- 
tween the two phases of large and small particles. This 
can be treated in an identical manner for polydisperse 
mixtures. 

This has been shown analytically for bidisperse sys- 
tems in [30] and validated experimentally in [33] . 

Generally, suppressing spurious numerical diffusion 
is a non-trivial task, requiring sophisticated finite dif- 
ferencing schemes [34] to maintain hyperbolicity. This 
type of stability, which generally controls the accuracy 
of the numerical results, is not an issue for solutions 
obtained using cellular automata [5 . 
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Fig. 15 A crushable flow with two mechanisms: segregation and comminution. Initially, the system is homogeneous, being a 
polydisperse sample with initial grainsize distribution defined by oti = —2. Four cases are considered with varying kb and k s . 
For each case, three plots are shown. Top: Evolution of the average grainsize at every height over time. Bottom left: Cumulative 
grainsize distribution at three points in the flow, corresponding to the crosses in the above plot. Top is red solid line, middle 
is green dashed line and bottom is black dash-dotted line. Bottom right: Cumulative grainsize distributions at steady state at 
same heights as previously. Solid blue lines represent best fit cumulative log-normal distributions. 



5.3 Segregation and comminution 

In many situations, segregation and comminution oc- 
cur simultaneously in a flow situation, such as in the 
grain milling depicted in Figure [I] In other cases, it 
is not even clear if segregation has occurred, yet par- 
ticles are advecting in space and strong comminution 
is observed, such as in earthquake faulting and snow 
avalanches. In many of these cases, we observe log- 
normal grainsize distributions, rather than power law 
distributions, which exist at all depths of flow. The ex- 
istence of these curves represents the competition be- 
tween the two mechanisms, where the comminution at- 
tempts to form a power law distribution, and the seg- 
regation attempts to create a locally monodisperse dis- 
tribtion. A log-normal distribution is one that obeys 
the following scaling for the cumulative grainsize dis- 
tribution Flat, 



"LN 



-erfc 
2 



Ins — /i 

aV2 



(19) 



where fi and a are the location and scale parameters, 
and erfc is the complimentary error function. 

In Figure [l5j simulations are shown in which both 
segregation and comminution are present. For all cases, 
log-normal cumulative grainsize distributions are ob- 
served over all depths, with p- values in the bottom half 
generally less than 0.001. 

As expected, increasing k s or kb decreases the time 
to reach a steady state in terms of the average grainsize 
s. It is evident that significant changes in the grainsize 
distribution will not occur indefinitely, even though seg- 
regation brings together particles of similar size, and 
comminution is accelerated by the segregation. 



5.4 Segregation, mixing and crushing 

We can now couple all three mechanisms and observe 
the evolution of the grainsize distribution as all of the 
constituent mechanisms interact. Each time step, we 
first check each cell and if the breakage rule is met, 
we change the cell's grainsize. Secondly, we iterate over 
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Fig. 16 A crushable flow with all three mechanisms. For all 
cases k s = kb = 1 and if = 1 is the height of the system. 
Initially, the system is homogeneous, being a polydisperse 
sample with initial grainsize distribution defined by ai = —2. 
Each row represents the same initial condition but with vary- 
ing amounts of mixing. From top to bottom, D = 0, 0.005 
and 0.05 respectively. Left: Cumulative grainsize distribution 
averaged at three heights z = 0.1, 0.5 and 0.9 at steady state 
(black, green and red respectively). Right: The average grain- 
size at every height evolving over time. 



all of the cells and swap them with a neighbour if the 
segregation rule is met. Finally, we again iterate over all 
cells and if the diffusion rule is met, we swap randomly 
with a neighbour. 

We now have a system that models avalanche and 
landslide flow, where particles at the base are sheared 
and crush, creating a lubrication layer, such as in Fig- 
ure 



15 The inclusion of mixing in the system, as shown 



in Figure [16j enhances the spread of sizes produced by 
comminution, and reduces the size of particles in the 
bottom-most layer of flow, enhancing the lubrication 
effect. Again, log-normal cumulative grainsize distri- 
butions are measured, which represent those found in 
many geophysical processes, such as in snow avalanches 



5.5 An equivalent continuum model 

Considering conservation of mass alone, we can express 
all three mechanisms in a continuum form as 



theories to represent internally (i.e. within the repre- 
sentative volume element) spatially correlated material. 
For the mechanisms of segregation and mixing, the 
length scale representing the local neighbourhood is not 
an important consideration. However for comminution, 
it must be included as part of the model to enable us 
to predict the correct final distribution. 



6 Conclusions 

We have shown that cellular automata can successfully 
capture the most important physical foundations be- 
hind evolving phenomena in crushable granular flows. 
To do this, we have used three distinct cellular au- 
tomata to explain the dominant mechanisms during 
such flows: comminution, segregation and mixing. 

By assembling all three cellular automata together 
we were able to explore the interactions between these 
phenomena. One surprising outcome is that in closed 
systems, crushable granular material are limited by power 
laws, however during flow the interaction with segrega- 
tion and remixing the system is limited by log-normal 
distributions. 

This paper highlights the power of the cellular au- 
tomata as a means to inspire continuum models. We 
have demonstrated that it is often possible to recover 
an analogous continuum description from the limit of 
the cellular automata rules. An interesting result of this 
tactic is the result that comminution does not follow 
this rule. For comminution, the ability to model cycles 
of crushing using a conventional continuum model with 
a grainsize coordinate is currently impossible, as some 
local rules are inherent to the system that must be in- 
cluded to describe the local configuration of grains. 

The success of this cellular automata is that it en- 
ables us to study the evolution and limits of the grain- 
size distribution in different scenarios. Current contin- 
uum models require a priori knowledge the grainsize 
distribution, and cannot educate us on the physical 
mechanisms involved in reaching this final state. 

IE acknowledges grant DP0986876 from the ARC. 
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This model differs from the cellular automata in the 
sense that the internal coordinate s does not retain 
information about local neighbours. Because of this, 
cycles of crushing cannot be produced. This omission 
gives an important insight into the use of continuum 
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